P2480 [SDOI2010]古代猪文

gdnCnn/dmod999911659g^{\sum_{d|n}C_{n}^{n/d}} \mod 999911659

gdnCndmod999911659g^{\sum_{d|n}C_{n}^{d}} \mod 999911659

因为 999911659999911659 是一个质数,由欧拉定理推论得

gdnCndmod999911658mod999911659g^{\sum_{d|n}C_{n}^{d} \mod 999911658} \mod 999911659

现在考虑求出指数部分,模数为 999911658=2×3×4679×35617999911658=2 \times 3 \times 4679 \times 35617

发现质因数的指数都为 11 , 其实就是特殊的扩展卢卡斯定理。只需要对每个质数分别计算后用中国剩余定理合并

#include <cstdio>
// #define int long long

const int MAXN = 35617 , Mod = 999911659;
int n , g;
int a[ 5 ] , p[ 5 ]= { 0 , 2 , 3 , 4679 , 35617 };

int Quick_pow( int x , int po , int p = Mod ) {
	int Ans = 1;
	for( ; po ; po >>= 1 , x = 1ll * x * x % p )
		if( po & 1 ) Ans = 1ll * Ans * x % p;
	return Ans;
}
int Exgcd( int a , int b , int &x , int &y ) {
	if( !b ) {
		x = 1 , y = 0;
		return a;
	}
	int r = Exgcd( b , a % b , y , x );
	y -= x * ( a / b );
	return r; 
}
int Inverse( int x , int p ) {
	int u , v;
	Exgcd( x , p , u , v );
	return ( u % p + p ) % p;
}
int Fac[ MAXN + 5 ] , Inv[ MAXN + 5 ];
void Init( int p = Mod ) {
	Fac[ 0 ] = 1;
	for( int i = 1 ; i <= p - 1 ; i ++ )
		Fac[ i ] = 1ll * Fac[ i - 1 ] * i % p;	
	Inv[ p ] = 0; Inv[ p - 1 ] = Inverse( Fac[ p - 1 ] , p );
	for( int i = p - 1 ; i >= 1 ; i -- )
		Inv[ i - 1 ] = 1ll * Inv[ i ] * i % p;
}
int C( int n , int m , int p = Mod ) {
	if( n < m ) return 0;
	return 1ll * Fac[ n ] * Inv[ m ] % p * Inv[ n - m ] % p;
}
int Lucas( int n , int m , int p = Mod ) {
	if( m == 0 || m == n ) return 1;
	return 1ll * C( n % p , m % p , p ) * Lucas( n / p , m / p , p ) % p;
}
int CRT( int *a , int *p ) {
	int Ans = 0 , Ps = 1 , P[ MAXN + 5 ];
	for( int i = 1 ; i <= 4 ; i ++ ) Ps *= p[ i ];
	for( int i = 1 ; i <= 4 ; i ++ ) {
		P[ i ] = Ps / p[ i ];
		Ans = ( Ans + (__int128)P[ i ] * Inverse( P[ i ] , p[ i ] ) % Ps * a[ i ] % Ps ) % Ps;
	}
	return ( Ans + Ps ) % Ps;
}

int main( ) {
	scanf("%d %d",&n,&g);
	if( g == Mod ) return printf("0") & 0;
	for( int i = 1 ; i <= 4 ; i ++ ) {
		Init( p[ i ] );
		for( int j = 1 ; j * j <= n ; j ++ )
			if( n % j == 0 )
				a[ i ] = ( a[ i ] + Lucas( n , j , p[ i ] ) + ( j * j == n ? 0 : Lucas( n , n / j , p[ i ] ) ) ) % p[ i ];
	}
	printf("%d\n", Quick_pow( g , CRT( a , p ) ) );
	return 0;
}